The Stability of Radiatively Cooled Jets in Three Dimensions 



Jianjun Xu^ 

Department of Astronomy, The University of Maryland, 
College Park, Maryland, MD 20742 

Philip E. Hardee 
Department of Physics & Astronomy, The University of Alabama, 
Tuscaloosa, AL 35487; hardee@athena.astr.ua.edu 

and 

James M. Stone 
Department of Astronomy, The University of Maryland, 
College Park, Maryland, MD 20742; jstone@astro.umd.edu 

Received March 2000 ; accepted June 2000 to appear ApJ, 543, November 10 



^IMS, MCIWorldcom, Vienna, VA 22182 



- 2 - 



ABSTRACT 



The effect of optically thin radiative cooling on the Kelvin-Helmholtz instability 
of three dimensional jets is investigated via linear stability theory and nonlinear 
hydro dynamical simulation. Two different cooling functions are considered: radiative 
cooling is found to have a significant effect on the stability of the jet in each case. 
The wavelengths and growth rates of unstable modes in the numerical simulations are 
found to be in good agreement with theoretical predictions. Disruption of the jet is 
found to be sensitive to the precessional frequency at the origin with lower frequencies 
leading to more rapid disruption. Strong nonlinear effects are observed as the result of 
the large number of normal modes in three dimensions which provide rich mode-mode 
interactions. These mode-mode interactions provide new mechanisms for the formation 
of knots in the flows. Significant structural features found in the numerical simulations 
appear similar to structures observed on protostellar jets. 



Subject headings: galaxies: jets — hydrodynamics — instabilities — ISM: jets and 
outflows 
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1. Introduction 



High resolution imaging surveys have revealed the association between Herbig-Haro (HH 
objects and protostellar jets ( [Reipurth et al. 1986 ; Strom et al. 1986; Mundt et al. 1987| ; Ray 



1987). Although individual HH objects can be formed via different mechanisms, most HH objects 



appear as bright knots along collimated jets, or associated with the working surface at the head 
of the jet, e.g., HH 4 6/47 ([Graham fc Elias 1983|), HH 4 ( |Mundt et al. 1984| ), HH 34 ( [Biihrke 



Mundt, k Ray 198j) , and HH 111 ( [Reipurth et al. 1997] ). Recent HST observations of HH 111 



( Reipurth et al. 1997| ) indicate that the bright knots along the jet beam arise from "bow" shocks 



associated with dense knots inside a collimated outflow. Sinusoidal distortion of the jet in 
HH 111 and staggering of the "bow" shocks suggest that the jet beam is helically twisted. HST 



observations of HH 30 ( Burrows et al. 1996 ) also reveal bright knots inside the jet beam. The 
proper motion of these knots has been detected and the observed bending of the jet changes with 
time, again suggesting helical twisting of the jet beam. 

A supersonic astrophysical jet is Kelvin-Helmholtz (K-H) unstable [see Birkinshaw (1991) for 
a review], and complex structures can be formed through the growth of instabilities in the jet 
beam. The study of the stability of supersonic astrophysical jets has focused on two fundamental 
questions: (1) What is required to make jets sufficiently stable to propagate large distances (up to 
a pc in the case of protostellar jets and hundreds of Kpc in extragalactic sources), and (2) what 
emission features, e.g., individual HH objects, are formed through instabilities in the jet beam? In 
protostellar jets, optically thin radiative cooling can strongly affect the dynamics, e.g., Blondin 
et al. (1990), Stone & Norman (1993a, 1993b, 1994), de Gouveia dal Pino & Benz (1993). These 
earliest numerical studies focused on propagation of the jets and interactions at the jet front, and 
showed that radiative cooling yielded strong fragmentation in a dense shell of material swept up 
by the working surface at the jet head. 

Interactions behind the jet front and the development of knots in a jet beam were first 
considered by Falle et al. (1987) who proposed that standing shocks could form the knots. 
Subsequently the knots were observed to have large proper motions (e.g. Eisloffel & Mundt 1995), 
and to have a close spacing inconsistent with standing shocks in the high Mach number flows. 
Later work showed that pulsation (temporal variation) of the outflow at the 10% level could 
produce aligned knots with bow shock morphology, e.g., Reipurth et al. (1986), Raga & Kofman 
(1992), Stone & Norman (1993b) de Gouveia Dal Pino & Benz (1994); Biro & Raga (1994); Biro 
(1996); Suttner et al. (1997). It was also suggested that knots can be the result of K-H flow 
instability, e.g., Biihrke, Mundt, & Ray (1988), Ray & Mundt (1993), Bodo et al. (1994), Rossi et 
al. (1997); Downes & Ray (1998). Additionally, the K-H instability could produce helically twisted 
structural features. That asymmetric structural features could be produced by K-H instability i s 
suggested by recent 2D linear sta bility analysis and 2D numerical simulations ( Hardee fc Ston^ 



1997 ; Stone, Xu, Sz Hardee 1997 ). This work showed that radiative cooling can increase the 
growth rate of K-H unstable modes significantly, and revealed a variety of asymmetric structural 
features formed as a result of instability. 

In this paper we study the K-H instability of 3D radiatively cooling jets through a systematic 
comparison between numerical simulations and linear stability theory over a wide parameter 
range. In §2 we describe the protostellar environment and radiative cooling assumed in this study. 
In §3 the 3D linear stability theory is summarized, and representative normal mode solutions 
appropriate to the numerical simulations are obtained. The numerical simulation results are 
shown in §4. In §4.2 we consider development of the surface wave modes and in §4.3 we consider 
development of the body wave modes and wave-wave interaction, and the pressure and velocity 
fluctuations arising from wave-wave interactions are compared to computations made from the 
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linear theory. Our numerical results are summarized in §5 and in §6 structures seen in the 
simulated jets are compared with observed protostellar jet structures. 



2. The Protostellar Jet Environment 

Protostellar jets show the same type of spectrum as low excitation HH objects, and the 
dominant emission lines include [S II] and the Balmer lines. The spectra are best explained 
as originating from a radiative shock, with velocity in the range 30 - 60 km s-i ( |Ray 1997] ) 



Although there are some YSO jets that contain high excitation emission, this emission can usually 
be identified with strong bow shocks where portions of the jet are moving into the ambient 
medium with a velocity of a few hundreds of km s~^, comparable to the typical velocities of 
~ 200 km s~^ exhibited by the fastest moving knots in the jet beam. The emission line spectra 
observed in typical YSO jets indicates that temperatures in the radiatively cooling regions are 
3 X 10^ — 10^ K, and line ratios, e.g., [S II] A = 6716A & 673lA, suggest that electron densities 
are Ng ~ 100 — 10, 000 cm~^ in the radiatively cooling regions. 

In addition to observational evidence that protostellar jets are overdense with respect to their 
surroundings (e.g., Ray 1997| ), results from numerical studies have shown that a propagating jet 



creates a cocoon behind the bow shock filled with hot gas that is less dense than th e jet. Dense 
molecular outflow s which are often observed to be associated with protostellar jets ( [Reipurth fc 



Cernicharo 1995| ) may surround the cocoon. The kinematic model of a dense jet embedded in a 



less dense cocoon which is in turn ensheathed by a very dense and cold molecular gas is consistent 
with detailed observation of a number of systems, e.g., HH 111 whose jet is surrounded by a 



bubble within the larger molecular outflow ( Nagar et al. 1997 ; Cernicharo fc Reipurth 199 



Since we focus our study on the stability properties of the jet beam which is only influenced by 
the jet's immediate surroundings, we can model a protostellar jet as a dense jet beam embedded 
in a hotter and less dense cocoon. In particular, we will assume that the protostellar jet is 
surrounded by a uniform, optically thin medium with low density (typically ~ 100 cm~^) and high 
temperature (typically ^ 10^ K). The fact that protostellar jets remain well collimated for lengths 
which are much larger than their radii argues they must be in pressure equilibrium with their 
surroundings. Numerical studies also have shown that the jet beam remains cold, neutral, and in 
approximate pressure balance with the cocoon medium. Pressure balance implies a jet temperature 
of Tjt = {nex/n'jt)Tex, where the subscripts ex and jt denote the externa l cocoon and jet gas 
respectively. Observations of some HH jets ( |Biihrke, Mundt, &: Ray 1988D indicate that typically 
the internal jet density exceeds the cocoon density by at least a factor of 10. For 10^ K and 

njt/nex = 10, Tjt ~ 10^ K. This indicates that the linear analysis of the equilibrium jet can focus 
on cooling processes which are effective with jet temperatures ~ 10^ K and cocoon temperatures 
^ 10'^ K. Because shocks with Vs ~ 100 km s~^ give temperatures as high as 10^ K, temperatures 

higher than 10^ K are considered in the cooling function in the numerical simulations. 

Optical emission from a protostellar jet and from the external gas results in a loss of internal 
energy from the system, and the loss of internal energy through radiative cooling can change the 
jet dynamics substantially. Thus, an important aspect of the numerical simulation of radiatively 
cooling jets is an accurate treatment of the microphysical heating and cooling rates. While it 
has been shown that a non-equilibrium ionization formalism significantly improves the accuracy 
of cooling terms (Stone & Norman 1993a), incorporating a time-dependent ionization fraction 
into a linear stability analysis is intractable. Thus, the cooling rate is assumed to be given by 
equilibrium cooling curves. In the work to follow two separate cooling curves are adopted: (1) the 
cooling curve for interstellar gas appropriate to protostellar jets calculated by Dalgarno & McCray 
(1972, hereafter DM), or (2) the curve described by MacDonald & Bailey (1981, hereafter MB) for 
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photoionized gas of reduced metallicity such as that found in elhptical galaxies and appropriate to 
extragalactic jets. While not strictly applicable to protostellar jets, the MB cooling curve allows 
us to explore the effect of temperature dependence of the cooling curve on jet dynamics. The full 
DM and MB cooling curves used in the simulations along with the piecewise power-law fits used 
in the linear analysis can be found in Figure 1 of Hardee & Stone (1997). 



3. The Linear Theory 

3.1. Normal Mode Analysis 

The stability of an adiabatic 3D jet with "top hat " profile residing in a uniform medium 
has been thoroughly investigated in the literature (see Birkinshaw 199 1] ). Modifications to the 
adiabatic 3D theory to account for radiative cooling are identical to those found in the 2D theory 
( Hardee fc Stone 19971) and we sketch the results here. The hydrodynamic equations of continuity, 
momentum, and energy are linearized within the jet and the external (cocoon) gas where u = 
and the flow velocity is reintroduced when solutions are matched at the jet boundary. The 
linearized hydro dynamical equations relevant to our model become 



^ + V.(po^'i) = 0, 



d_ 
dt 
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where, in general, the perturbed quantities are written as p = po + pi, v = u- 
The energy equation is written following Hunter & Whitaker (1989) with 
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where Cq is the initial cooling rate, a = (FPo/po)^^^ is the sound speed, and T is the adiabatic 
index. 



In cylindrical geometry a random perturbation of pi, vi, and Pi to an initial equilibrium 
state Po, u, and Po can be considered to consist of Fourier components of the form 



/i(r, (p, z) = fi{r) exp[i(A;2; ± n(/) - ujt)] 



(4) 



where flow is along the z-axis, and r is in the radial direction with the fiow bounded hy r = R. In 
cylindrical geometry k is the longitudinal wavenumber, n is an integer azimuthal wavenumber, for 
n > the wavefronts are at an angle to the fiow direction, the angle of the wavevector relative to 
the flow direction is 6 = tan(n/kR), and +n and —n refer to wave propagation in the clockwise and 
counterclockwise sense, respectively, when viewed outwards along the fiow direction. In equation 
(4) n = 0, 1, 2, 3, 4, etc. correspond to pinching, helical, elliptical, triangular, rectangular, etc. 
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normal mode distortions of the jet, respectively. Propagation and growth or damping of the 
Fourier components is described by a dispersion relation 

Xjt JniPjtRjt) Xex HniPexRjt) 

where the primes denote derivatives of the Bessel (J) and Hankel (H) functions with respect to 
their arguments, /3 and x are given by 



Pex — [ + 9 Qez] ^ 1 



(ijt = {-k H ^2 ■ , 

and 



_ 2 

Xex — Pex^ ) 



^2 



Xjt = Pjt{uJ - kuY 

All heating and cooling information is contained in the dimensionless terms Qjt and Qex where 
and 

_ 1 + iF^y{uj - ku) 
~ 1 + iF^'/{uj - ku) ■ 

For normal mode n the axial wavelength associated with a 360° helical twist of a wavefront around 
the jet beam is given by = n\n where = 271 /k. The angular frequency appearing in the 
linear analysis, w, represents a 360°/n helical twist and is properly related to an angular precession 
frequency by a; = nujp. For example, if an elliptical jet distortion, normal mode re = 2, rotates 360° 
with the precession frequency oOp, at a fixed azimuthal angle the frequency at which the jet surface 
oscillates is 2ojp. In general, each normal mode, re, consists of a single "surface" wave and multiple 
"body" wave solutions to the dispersion relation. Structural differences between the surface and 
body waves are discussed in §3.4. 



3.2. Heating and Cooling Rates 



Following previous 2D work ( [Hardee &: Stone 1997 ; ptone, Xu, &: Hardee 1997| ), we chose a 



jet of radius Rjt = 2.5 x 10^^ cm with number density rijt = 600 cm~'^, temperature Tjt = 10^ 
K and sound speed ajt = 3.73 x 10^ cm s~^. Initially the jet is in pressure equilibrium with an 
external (cocoon) gas with number density Uex = 60 cm~^, temperature Tg^ = 10^ K and sound 
speed ttex = 1-18 x 10^ cm s~^. With heating and cooling rates of the form H = Aun^^T^" , and 
C = Acn"^T^^ , Fp and Fp are given by 

Fp= ^^^~^^^ C{f3c-f3H), 
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and 

Fp = -^C [{Pc - Ph) + {an - ac)] , 

where C = Acn'^^T^^ ergs cm-^ s"^. We will assume that the heating rate is independent of 
the temperature, i.e., Ph = 0, and is proportional to the density, i.e., an = 1, and the initial 
equilibrium requires that Hq = Cq. 

In the linear analysis a piecewise power-law fit to the DM cooling function is represented by 

1 3.0 X 10-28T0-5 erg cm^ s'^ : T < 10"^ K 

[ 1.5 X 10-242^0.5 s~i . 2^ > iqA 

For DM coohng, C = Ac-n^^Tg^^ = uIAdm and for our choice of parameters 

5 . „2r0.5 _ / 2.16 X 10-9 s-1 \ 




and 




9poa2 ^ ^ I 1.37 X 10-11 s-i 



A „2^o.5 _ j -1-29 X 10-9 s-i 



3poa^ [-0.82 X 10-11 s-ij' 

A power-law fit to MB cooling curve of the form 

A^^^ = 1.5 X 10-342^2.53 gj.g ^^3 g-l 

is used in the linear analysis to simulate the MB cooling curve. For MB cooling, 
C = Ac'tiq'^Tq'^ = HqAmb and for our choice of parameters 

f-Fp""! 25 . . , „ f 1.46 X 10-10 s-il 



and 



Fn ^ 9poa2 ^ ° ° " \4.30 x lO"" s'l 



1 A 2^2 53 _ j 5-24 X 10-11 S-1 

^— AcnoJo -i 1.54 X 10-11 s-i 




The most important difference in jet stability properties between DM cooling and MB cooling 
arises from the change in sign of Fp ([Hardee &: Stone 1997]) . The positive value of Fp for MB 
cooling is a consequence of the steeper temperature dependence of the MB cooling function. 
Positive values of Fp occur when {f)c — Ph) + {oh — etc) > 0. For our choice of heating rate 
depending linearly on the density and of cooling rate depending on the density squared, Fp > 
when > 1- The DM power-law fit used in the linear analysis serves to illustrate the effect of 
a shallow dependence of radiative cooling on temperature in both jet and external fluid. The 
MB power-law fit used in the linear theory serves to illustrate the effect of a steep dependence 
of radiative cooling on temperature in both jet and external fluid. In the simulations and in the 
theory the "equilibrium" heating rate Hq = Cq is determined by the requirement that the jet and 
external medium be in thermal equilibrium initially. Because cooling rates are different in the 
jet and external fluids as a result of temperature and density differences, the initial heating rate 
required to establish and maintain thermal equilibrium is different in the two fluids. 



- 8 - 



3.3. Dispersion Relation Solutions 

Numerical solutions to the dispersion relation for adiabatic (AD), Dalgarno &; McCray (DM) 
cooling and MacDonald & Bailey (MB) cooling jets have been obtained for the pinch, helical, 
elliptical, triangular and rectangular normal modes for M^x = uja^x = 5 and 20 jets, jet speeds 
of 59 km s^^ and 236 km s~^, respectively. The solutions for M^x = u/ciex = 5 and 20 are 
qualitatively similar, and in Figure 1 we show those appropriate to Mex = 20 and a jet speed of 
236 km s~^. In the numerical simulations, and possibly in protostellar systems, jets are perturbed 
at their origin by a periodic motion at some frequency; thus, we solve the dispersion relation for a 
complex wavenumber as a function of angular frequency. A negative value for the imaginary part 
of the wavenumber indicates spatial growth at that angular frequency with an e-folding length 
~ \kj I • 

We identify the following features in the 3D solutions: 

(1) The solutions for the pinch and helical normal modes in 3D are very similar to the solutions 
obt ained in 2D for the sym metric (pinch) and asymmetric (sinusoidal) normal modes of the slab 
jet ( [Hardee fc Stone 19971) . The higher order normal modes of the 3D jet which have no analog 
on the 2D jet, e.g., the elliptical etc. normal modes, can be thought of as harmonics of the helical 
mode and behave similarly to helical mode solutions. 

(2) In general, the linear growth rates, scale inversely with the Mach number (not shown), 
and the linear growth rates, especially for the body modes, are relatively smaller in 3D than in 
2D. Body wave mode growth rates can be larger, comparable to or smaller than the growth rate of 
the surface wave mode for the pinch, helical and elliptical normal modes, respectively. Body wave 
mode growth rates are less than the surface wave modes for all higher order normal modes (not 
shown) . 

(3) The effect of radiative cooling on the dispersion relation is similar in 2D and 3D, i.e., DM 
cooling significantly increases the linear growth rate at higher frequencies and MB cooling 
decreases the linear growth rate at higher frequencies. The high frequency growth rate plateau 
for DM cooling extends to frequencies only slightly higher than those shown in the figure, at 
which point the growth rate rapidly declines towards the adiabatic results at high frequency. 
Adiabatic and MB cooling jets show distinct "resonant" (fastest growing) angular frequencies, 
Lv*, and wavelengths, A*, for almost all surface and body wave solutions, and the DM cooling 
jet shows distinct resonances for the body wave solutions. As was found in 2D the presence of 
DM cooling adds a pinch cooling mode (Ps2) that does not exist for AD or MB cooling jets. In 
general, the wave modes are purely real on the adiabatic jet or damped on the radiatively cooling 
jets (damping rates not shown) in regions where they are not growing. 

3.4. Fluid Displacements, Velocity & Pressure Fluctuations 

Displacements, $,{r,(j),z), of jet fluid from an initial position {r,(p,z) can be written in the 

form 

^(ro, zs) = A(ro)e*^(^«)er(i?) exp[f(A:z, ± - tot)] , (6) 

where Zs and (ps are the axial and azimuthal positions at the jet surface, ro is the initial radial 
position, S,r{R) is the radial displacement at the jet surface, the A(r)e*^^^^ are given by equations 
(AlO) in Hardee, Clarke, &: Rosen (1997, hereafter HCR), and the k{uj) are normal mode solutions 
to the dispersion relation. Fluid displacements are modified in amplitude and rotated in azimuthal 
angle or shifted along the jet axis relative to those at the jet surface by A(r)e*'^^^\ The difference 
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in structure between surface and body wave solutions on adiabatic jets at their resonant frequencies 
as revealed by displacement surfaces is shown explicitly by Figure 13 in HCR. The accompanying 
velocity perturbation, vi(r, 0, z) = d$,/dt, and pressure perturbation, Pi{r,(j),z), can be written in 
the form 

vi(r', (/)', z') = -i{LO - ku)A{ro)e'^'^'^'>^Cr{R) exp[i{kzs ± n</., - ut)] , (7) 
Pi(r', (/.', z') = B{r^)e'^^^'-^^UR) exp[i(A;z, ± - ujt)] , (8) 

where points on initially cylindrical surfaces are displaced radially to r' = vq + Sr, axially to 
z' = Zs + 6z, and azimuthally to 4>' = (ps + 6(1), and 6r, 6z, 6(p are the components of ${ro, (ps, Zg)- 
In equation (8) B{ro)e'^p(-''-<^^ = {xjt/(3jt)[Jn{Pjtro)/JUPjtR)] [e.g., Pardee et al. 199^ eq. (14)]. 



The total velocity is given by v = u + vi and the total pressure is given by P = Pq + Pi. 

If the azimuthal and axial phase shift is small then the radial fluid displacement of a "surface" 



wave mode n > in the jet interior at constant azimuthal angle is Cr{^) ~ ^r{R){f / R)^ ^ ( Hardee 



1983 ). The accompanying velocity and pressure variations produced by higher order surface modes 
also show a rapid decrease inwards. On the other hand, at a constant azimuth the "body" wave 
modes have a reversal in fluid displacement at "null displacement" surfaces in the jet interior and 
typically the maximum pressure is near the null surface. 

Jet distortion and accompanying pressure and velocity structure associated with individual 
normal mode surface and body waves at frequencies les s than or equal to t he resonant fr equency 



have been thoroughly investigated for adiabatic flows ( [Hardee et al. 199S ; Hardee 20001) . Jet 



distortion and the accompanying pressure and velocity structure across the jet and parallel to 
the jet axis are not strongly modified by radiative cooling and the structure of the individual 
wave modes on adiabatic and cooling jets as a function of frequency at and below the resonant 
frequency of the appropriate mode is similar. The numerical simulations presented in the next 
several sections perturb the jets at low frequencies and at high frequencies above the resonant 
frequency. At the higher frequencies the simulations reveal interesting wave-wave interactions that 
have not been investigated previously. We will consider the high frequency structure of various 
modes in comparison with structures observed in the numerical simulations in §4.3. 



4. Numerical Simulations 



4.1. Initial and Boundary Conditions 

All of the simulations were performed using the modified three-dimensional hydrodynamic 
code CMHOG. For details on the hydrodynamical algorithms see Stone, Xu, & Hardee (1997). 
At dimensionless simulation time Tgim = {O'ex/Rjtj't = 0, a jet beam of Mach number M^x with 
uniform density rijt and z-velocity = M^xO^ex is set up across a Cartesian computational 
grid in the z direction. Outflow boundary conditions are used except where the jet enters the 
computational grid where inflow boundary conditions are used. The size of the computational 
domain is varied in the axial direction depending on the Mach number of the jet and the 
perturbation frequency in order to ensure that many wavelengths of the most unstable mode are 
captured. Simulations with DM and MB type cooling were performed, along with simulations of 
adiabatic jets (AD) that serve as a baseline for comparison. Grid size along with all the other key 
parameters associated with the simulations are listed in Table 1. Note that the duration of various 
simulations, Tsim, is different. 

In these simulations we study the dynamical evolution of unstable jets that are initially 
in a delicate pressure equilibrium with a low-density ambient medium and in which the gas is 
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established in thermal equilibrium with net cooling rate n^A — nH = 0. To maintain thermal 
equilibrium a different heating rate, H, is used in jet and ambient fluid and the initial equilibrium 
could be maintained for many dynamical times if the system was not perturbed. As the perturbed 
simulations evolve and the fluids mix the heating rate, H, is varied according to whether a grid 
zone temperature is above or below 10,000 K for DM cooling or according to the fraction of ambient 
and jet material in a grid zone for MB cooling. These different approaches reflect the discontinuity 
in the DM cooling function at 10,000 K and the continuity of the MB cooling function, and are 
identical to the methods used in the previous 2D simulations. The per particle cooling rate. A, is 
computed in each grid zone according to the complete DM or MB cooling curves (see Figure 1 in 
Stone, Xu, & Hardee 1997). Note that the net cooling rate is now non-zero but is less than nPA. 



As a result the cooling length (cool > Vshocktcooi ~ Vjttcooi where tcooi = (r - 1) ^{kT/nA), and 
shocked regions in the ambient medium can be underesolved with cooling length a few times the 
grid spacing. However, underesolution of these shocks does not influence the stability or internal 
structure of the jet beam which is primarily what we study here. 

A periodic precession of the jet velocity is applied at the jet nozzle to break the symmetry. In 
the simulations, the pitch angle, 0, of the jet velocity relative to the axial direction varies from 9 



0.0025 to 0.01 rad d epending on the precession frequency. As in previous 2D simulations ( Stone 
Xu, fc Hardee 1997 ) we choose the initial precession (perturbation) frequency to be ojRjt/u = 0.1 
0.4 and 1.0 for Mach 5 jets, and ujRjt/u = 0.025 and 1.0 for Mach 20 jets. With our assumed jet 
radius and velocities the precession period ranges from about 840 to 21 years. 



4.2. Surface Wave Modes 



4-2.1. Low Frequency Precession 

We performed simulations with angular frequency ujRjt/u = 0.1 and LoRjt/u = 0.025 on 
Mach 5 and 20 jets, respectively. This angular frequency is about a factor of three below the 
resonant frequency on the AD jet, and for our assumed jet radius and flow speeds corresponds 
to a precession period Tp = 843.7 yr. Mach 20 jets have speeds that are comparable to observed 
protostellar jets. As comparable behavior is seen in the Mach 5 and 20 numerical simulations, 
albeit with different length scale we show a volumetric rendering of the density for Mach 20 AD 
and DM jets in Figure 2. At precession frequencies much below the AD helical resonance, the 
linear theory indicates that helical and higher order surface wave modes will be dominant as the 
surface mode growth rates are much higher than the body mode growth rates on AD jets and 
body modes are damped on DM and MB jets at this low precession frequency (see Figure 1). 

The dominant oscillation is the result of helical twisting excited by the precession at the 
inlet. Note that the Mach 20 AD jet remains collimated across most of the computational grid 
even though the linear growth rates of AD and DM jets are nearly the same at this low precession 
frequency. Thus, non-linear processes speed the breakup of radiatively cooling jets relative to 



adiabatic jets. A si milar result was found in 2D for identical parameters (see Figure 9 in Stone 
Xu, fc Hardee 1997|) . Some high density knot formation is apparent in the region where the Mach 



20 DM jet breaks up at a distance slightly less than 400i?j(, approximately 2.3 e-folding lengths. 
A 2D DM jet breaks up at a slightly longer distance at this precession frequency even though the 
2D sinusoidal surface mode growth rate is larger than the 3D helical surface mode growth rate 
(about 50% larger). Upon break up the 2D jet develops denser knots. We note that the Mach 5 
DM jet develops similarly but only propagates about 1/4 the distance of the Mach 20 jet ~ about 
lOOiJjt at break up. This difference is about the factor indicated by the direct scaling between 
growth length and Mach number predicted from theoretically computed spatial growth rates. Very 



- 11 - 



little knot formation occurs at the lower speed of the Mach 5 jet in 2D or in 3D simulations. 

In the region where Mach 20 and Mach 5 DM jets break up the structure is very complicated 
and bow shocks are the dominant feature. This region also shows rapid mixing of jet and ambient 
material. The mixing, as measured by an entrainment volume and by an entrained mass, for the 
Mach 5 DM simulation as a function of position along the jet beam at time Tsim = 100 is shown in 
Figure 3. The entrainment volume is defined as that volume of the fluid that has C > e where 
< e < 1 (e.g., [Loken et al. 199q ), and C is a "color" variable that traces the jet fluid. A value 

= 1 indicates that all zones on the grid in a plane transverse to the jet axis at location z contain 
jet material with C > e. The entrained mass Ms is defined as the mass of fiuid with an axial 
velocity Vz > Suex (e.g., [Loken et al. 199(^ ; Bassett &: Woodward 199"^ ). Mixing is considerably 
enhanced in 3D relative to a similar 2D Mach 5 numerical simulation (see Figure 11 in [Stone, Xu, 
fc Hardee 1997 ). The volume containing entrained material saturates near z = 50Rjt in 3D versus 



SOORjt in 2D. This is due to the dramatic increase in area of the jet ambient medium interface 
when the 3D jet disrupts, thus accelerating the mixing process. The entrained mass shows an 
increase only for axial velocities of jet and entrained material that are less than the external sound 
speed, i.e., 6 < 1, whereas in 2D the entrained mass showed an increase for 6 < 2. This is a result 
of the fast break-up and rapid slow down of jet material accompanying the mixing process in 3D. 
Note that the decline in entrained volume and mass at axial distances z > 125Rjt has occurred as 
fiow from the inlet no longer reaches these distances at time r^jm = 100. Presumably the slowly 
moving mixed material would progress across the computational grid at longer times. 



4-2.2. Moderate Frequency Precession 

In the Mach 5 simulations the angular precession frequency ujRjt/u = 0.4 is on the order of 
the resonant frequency of the first few normal mode "surface" waves on the AD and MB jets, and 
is within the unstable frequency range for the important accompanying "body" waves on AD, 
DM, and MB jets. For our assumed jet radius and fiow speed this angular frequency corresponds 
to a precession period of Tp = 210.9 yr. No comparable simulations were performed for Mach 20 
jets. The jets showed helical twisting at wavelengths, X/Rjt ~ 14.3 (AD), 15.6 (DM), and 14.3 
(MB), very close to those predicted theoretically - 14.1 (AD), 15.6 (DM) and 14.0 (MB) - where 
theoretical wavelengths are given by X/Rjt = 2TTVgp/u! and Vgp = {du! / dk)\j.eai is the group velocity 
at the precession frequency. Thus, the linear analysis correctly predicts the helical wavelengths. 

At this precession frequency the surface modes all have substantial growth rates. Ultimately, 
growth of these modes can cause the jet to bifurcate or trifurcate. Splitting the jet into filaments 
provides a large interface between the jet and ambient medium, and the observed rapid mixing 
with the external medium. Quantitative comparison between linear growth rates shown in Table 2 
indicates that the linear growth rate for the helical surface mode is smallest for the adiabatic jet, 
larger for the MB cooling jet, and largest for the DM cooling jet. At least qualitatively, the more 
rapid predicted growth of helical instability on the cooling jets is verified by the simulations. The 
rapid growth of higher order surface modes (see §4.2.4) in 3D precludes a quantitative analysis of 
amplitude growth of the helical mode. In any event, at a moderate precession frequency, radiative 
cooling has accelerated the breakup of the jet beam, just as we found at a low precession frequency. 

Outside the jet beam the simulations show a pattern of spiral shocks produced by the 
large amplitude helical twist. In the external gas this shock, when viewed edge on, appears as 
alternate bow-shaped arcs on each side of the jet beam. This is illustrated in Figure 4 by a 
volumetric rendering of the temperature and of the density from the Mach 5 DM simulation. 
High temperatures lie immediately behind the shock front, and the volumetric rendering of the 
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density shows that the shocks originate at dense knots on the jet surface. In Figure 5, we show a 
volumetric rendering of p^T (density squared times temperature) for this Mach 5 DM cooling jet. 
The p^T image emphasizes regions in the jet beam where the jet material is shocked and raised to 
higher temperatures and simulates emission from a single line, e.g., [S II]. At this frequency, the 
helical ridge on the surface of the jet beam is caused by the helical surface wave. The "hot spots" 
are caused by projection combined with "limb brightening". The high temperature postshock 
gas behind the spiral shock forms a thin continuous radiating layer. If observed from a direction 
which is perpendicular to the jet beam, the radiative region looks like discrete "knots". At larger 
distances the jet starts to entrain ambient material and dense knots are formed via cooling. These 
knots have a velocity slightly less than the initial jet velocity and a temperature which is much 
higher than the initial jet temperature. 



4-2.3. High Frequency Precession 

The precession frequency ujRjt/u = 1.0, about three times the helical resonant frequency on 
the Mach 5 adiabatic jet, corresponds to a precession period of Tp = 84.4 yr for our assumed jet 
radius and flow speed. Figure 6 shows the jet density in a slice plane along the jet axis for AD, 
DM, and MB Mach 5 jets at this frequency. No comparable simulations were performed for Mach 
20 jets. The jet beams do not show a large amplitude helical twist, but the internal patterns 
reveal the presence of periodic knots inside the AD and DM jet beam. This result is similar to 



numerical simulations in 2D slab geometry ( Norman fc Hardee 1988 ; Hardee, Cooper, fc Clarke 
1994 ; ^tone, Xu, &: Hardee 1997| ) which have shown that the jet beam cannot respond bodily to 



perturbations much higher than the resonant frequency, i.e., the entire jet beam cannot form high 
frequency sinusoidal or helical patterns with significant surface displacement. 

While these jets do not exhibit a large scale helical twist associated with the surface wave 
at the precession frequency, the presence of a resonant wavelength in the AD and MB jets is 
exhibited by a helical twist seen in the AD and MB simulations as the jets break up. In particular, 
we measure a helical twist in AD and MB simulations with wavelength A ~ 15Rjt ± 2Rjt. The 
linear theory predicts a resonant wavelength for the helical surface wave of A* = 16.9Rjt and 
13.2Rjt for AD and MB jets, respectively, in good agreement with the simulations. Thus, jet 
breakup is associated with helical twisting at about the resonant wavelength. These jets remain 
well collimated to distances in excess of double that associated with the moderate precession 
frequency that was near to the helical surface wave resonance. Thus, the collimation of the jets 
is preserved fairly well due to the inability of jet beams to respond bodily to the high precession 
frequency and very small perturbation at the lower resonant frequency. 



4.2.4. Higher Order Surface Modes 

Even though the precession that we use preferentially excites the helical surface mode, the 
higher order surface modes have higher maximum growth rates on AD, DM, and MB jets (see 
Figure 1) and can grow to significant amplitudes. Growth of the higher order modes is observed to 
lead to elliptical, triangular, rectangular, etc. distortions of the jet cross section in AD, DM, and 
MB Mach 5 jet simulations. This type of distortion is seen in all Mach 5 simulations performed 
with precession frequencies ujRjt/u = 0.4 & 1.0. The lower of these two precession frequencies is at 
about the calculated resonant frequency for the helical through rectangular surface modes on AD 
and MB Mach 5 jets, is slightly below the high frequency growth rate plateau on the DM jet, and 
should excite these surface modes at about their maximum growth rates. The higher precession 
frequency lies in the range where body modes have their maximum growth rates and where the 
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surface modes have smaller growth rates and (or for the DM jet with high frequency surface mode 
growth rate plateau) have a lesser effect on jet distortion. 

The development of cross section distortion from the Mach 5 DM simulation with precession 
frequency tuRjt/u = 0.4 is shown in Figure 7. Qualitatively similar results are found for the 
comparable Mach 5 AD and MB simulations. As first observed by Hardee & Clarke (1995), we 
find that higher order modes with the fastest growth rates appear first, and as one moves down 
the jet beam away from the nozzle, the cross section area of the jet first becomes pentagonal, then 
rectangular, triangular, and finally elliptical. Note the jet beam bifurcation as a result of the 
growth of the elliptical mode. Quantitatively, the distance where we expect to see the higher order 
modes to appear will be proportional to £* = \kj where kj is the maximum spatial growth rate. 
With our small amplitude initial perturbation, growth to a significant amplitude should require a 
distance of about 27r/A;J, i.e., about six e- folding lengths, and these distances are listed in Table 2. 
We find that rectangular, triangular, and elliptical distortions appear in the same order along AD, 
DM and MB jets, and the locations of their appearance are close to the values listed in Table 2. 
Thus, the linear analysis correctly predicts the relative growth rates of the surface modes. 

At a precession frequency of ujRjt/u = 1.0, the Mach 5 AD, DM and MB jets remain well 
collimated to about twice the distance of their uiRjt/u = 0.4 brethren. Cross sections (not shown) 
reveal lesser surface distortion but with internal structure consistent with the development of 
body modes. Thus, high frequency precession, where the linear analysis predicts that surface 
mode growth is reduced and/or perturbs only a shallow surface layer of the jet beam, is shown 
by the numerical simulations to be much less destabilizing to the jet beam than lower frequency 
precession as surface modes are suppressed. 

4.3. The Body Modes & Wave- Wave Interactions 

4-3.1. Low Frequency Precession 

At precession frequencies much below the AD helical resonance, the linear theory indicates 
that helical and higher order surface wave modes will be dominant, and the body wave modes are 
either purely real on the AD jet or arc damped on the DM and MB cooling jets. While growth 
may occur for all modes near resonance, the absence of initial perturbations at higher frequencies 
suggests that growth of the surface and/or body wave modes at resonance will be slow. In fact we 
find evidence for interaction between the helical surface wave and a weakly damped helical body 
wave where both are excited by the processional perturbation. 

An interaction between helical surface and body waves is shown in Figure 8, which contains 
composite volumetric renderings of p^T from the Mach 5 DM cooling jet perturbed at a frequency 
of ujRjt/u = 0.1. The top panel shows the entire jet beam. The middle panel shows "hot spots" 
in the jet beam by excluding jet material with temperatures at or below 1, 000 K. The bottom 
panel shows an enlargement of the knot region. The hot spots have temperatures ranging from 
1000 K — 10,000 K and appear very prominently at locations of z ~ 23 & 37 Rjf However, at 
the precessional frequency the helical and higher order surface modes, have a wavelength of about 
58Rjt (from the linear analysis). Note the long wavelength oscillation evident in the images. On 
the other hand, the helical body modes are only weakly damped at this frequency with damping 
length ~ lOOORjf. If the helical surface wave with wavelength Xgurf = 58.0Rjt, and the first helical 
body wave with wavelength Xbody = 18.47Rjt are in phase at the origin, subsequent in phase 
locations should occur at 23Rjt and at 35.7Rjt. These positions are very close to the observed 
positions of the prominent knots at z 23 &; 37 Rjf. The interaction can also be seen in jet cross 
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sections (not shown). Thus, the interaction between body waves and surface waves can lead to 
knot formation. Knots or hot spots created in this fashion are stationary. 

A second type of interaction also appears possible. In the Mach 20 DM jet perturbed with a 
precessional frequency oi ujRjt/u = 0.025 the jet ambient medium interface appears to break up at 
several locations (see Figure 2 and compare the Mach 20 AD and DM jets in the first half of the 
computational grid). The break up of the jet-ambient medium interface in this simulation appears 
to be the result of energy deposition at the jet surface from a damped internal body wave. With 
the exception of this feature the Mach 20 DM jet shows very little structure or evidence for dense 
knot formation inside the jet beam. 

4-3.2. Moderate to High Frequency Precession 

The angular frequencies ujRjt/u = 0.4 & 1.0 are within the unstable frequency range for the 
important "body" waves on AD, DM, and MB Mach 5 jets. In the DM and MB simulations at 
the lower frequency, we observe dense filaments aligned with jet flow inside the jet beam near to 
the surface, and see indications of dense knot formation. The location and length of these dense 
filaments is coincident with constructive wave-wave interactions between the first three helical 
body waves if all waves are launched from the origin in phase. For example, density cuts in a slice 
plane along the jet axis from the DM and MB simulations (not shown) reveal a dense filament of 
length ~ A.lRjt at z < lORjt- We note that full wavelengths of the first, second and third helical 
body waves are \*/Rjt = 10.1, 7.5 and 6.0, and the first helical body wave at 2; = WRjt has 
over-run the third helical body wave by a distance of ~ i.lRjt with less over-run for the second 
helical body wave. This suggests that compression over this length interior to z = lORjt by the 
three body modes has led to radiative cooling and the observed filament formation. 

At the higher angular frequency of uRjt/u = 1.0 density cuts in a slice plane along the jet 
axis shown in Figure 6, and plots of pressure, axial and transverse velocities along the jet axis, 
shown in Figure 9 for the DM jet, reveal knot structure in the Mach 5 AD and DM jets with 
periodic spacing of 34.0Rjt and 27. QRjt, respectively. Similar structure is not seen on the MB jet. 
While the structure suggests periodic pinching and knot formation, the periodic spacing is much 
longer than the wavelength of any normal mode body or surface wave at the precession frequency. 
However, close examination of "resonant" wavenumbers for first and second helical body modes 
(Table 3) reveals a potential beat pattern with spacing of 39.2i?jt and 26.2Rjt for AD and DM 
jets, respectively, in excellent agreement with the observed knot spacing. The lack of a beat 
pattern and knots in the MB jet might be due to the low growth rate of the second body wave 
relative to the first body wave on the MB jet. 

How can periodic knots be produced by a beat pattern between helical body modes? Since 
only the beat pattern between the two helical body modes can generate the observed wavelength, 
we interpret the simulation results to mean that the wave-wave interaction between the helical 
body waves has excited and funneled energy into a pinch mode at the long beat wavelength. To 
show how this could be so we have computed the pressure and velocity fluctuations accompanying 
the first two helical body modes, Hbl and Hb2 using solutions computed from the dispersion 
relation for the Mach 5 DM jet at the precession frequency loR/u = 1, the fluctuations resulting 
from their beat pattern (Hbl2), and the fluctuations associated with the pinch surface cooling 
mode Ps2 at the beat wavelength. Fluctuations were computed using equations 6 - 8 in §3.4. 
The results for pressure fluctTiations appropriate to those observed in the Mach 5 DM simulation 
(Figure 9) are shown in Figure 10. In Figure 10 we plot component structure along ID cuts 
parallel to the jet axis (z-axis) at different locations on the transverse +?/-axis where velocity 
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components and Vy represent azimuthal and radial velocities, respectively. Coarse mapping of 
values computed in cylindrical {r,4>,z) coordinates to the Cartesian (x,y,z) coordinate location of 
the ID cuts has led to the raggedness of the lines in the figure. 

The ID cuts are shown at multiple radial locations between jet center (dotted lines) and 
surface (dashed lines) for the first (Hbl) and second (Hb2) helical body waves. The highest 
pressures associated with Hbl at this frequency are near the jet surface. The Vy (radial) velocity 
component shows opposing motions between material near the jet surface and in the jet interior, 
and a half-wavelength axial phase shift between outwards radial motion near the jet surface 
and at the jet center. The highest pressures associated with Hb2 are in the jet interior with a 
half-wavelength axial phase shift in pressure extrema between the jet surface and the jet center. 
At this frequency the radial velocity shows a large axial phase shift between jet surface and jet 
center. The axial velocity fluctuation is small for both helical body modes near jet center. The 
combination of the two helical body waves (Hbl2) leads to a combined fluid displacement different 
from that of the individual waves, and the resulting pressure and velocity fluctuations along the jet 
axis (dotted lines) for the individual waves cannot simply be added to give the pressure fluctuation 
(dash-dotted line) and beat pattern in the transverse velocity components along the jet axis. Still 
the difference between the wavenumbers of the first and second helical body modes (see Table 3) 
leads to the beat pattern with wavelength of ~ 2QRjt. The transverse velocity components along 
the axis (dashed line & dotted line) are 90° out of phase as expected for helical distortion. Note 
the lack of fluctuation in the axial velocity component as is observed in the simulation. Pressure 
and velocity fluctuation structure resulting from this wave-wave interaction changes significantly 
off the jet axis. 

The pressure and velocity fluctuations expected to accompany the pinch cooling mode (Ps2) 
between the jet axis and surface (dash-dotted line) at a wavelength comparable to the beat 
pattern between the two helical body mode waves and with a comparable pressure fluctuation 
shows a substantial axial velocity fluctuation and a lack of radial velocity fluctuation along the jet 
axis (dotted line) as is observed in the simulation. Pressure and axial velocity fluctuation do not 
change as a function of the radius but radial velocity fluctuation increases slightly towards the 
jet surface (dashed line). Clearly wave- wave interaction between the helical body modes couples 
via the pressure fluctuation near jet center to the pinch cooling mode with pressure fluctuation 
occupying a large fraction of the jet interior and with some non-linear interaction between the 
velocity components that results in both small axial and transverse velocity fluctuation near jet 
center. The resulting nearly axisymmetric periodic pressure and density fluctuation produces the 
observed knots. 



4-3.3. Very High Frequency Precession 

In the Mach 20 DM simulation we study jet response to an angular precession frequency 
of uRjt/u = 1.0 that is over ten times the AD jet helical resonant frequency. For our assumed 
jet radius and flow speed this angular frequency corresponds to a precession period of Tp = 21.1 
yr. At this frequency, the linear theory indicates a high frequency plateaTi in the surface wave 
mode growth rate that is above the body wave mode growth rates (see Figure 1). The helical 
nature of the initial precession guarantees that modes higher than helical are unlikely to develop 
to significant amplitude, although growth rates are comparable. We note that the wavelengths 
predicted to accompany this precession frequency are shorter for the helical surface mode than for 
helical body modes. This is a reversal of the usual relationship. 

This simulation provides another interesting example of knot formation via wave-wave 
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interaction. Knot formation appears in the p^T volumetric rendering image in the upper panel 
in Figure 11 at Tgim = 5. The volumetric rendering reveals a clear helical pattern within the 
jet beam with wavelength 6.5Rjt- From the linear analysis, we calculate wavelengths of Xh = 
A.SRjt, 6.65Rjt, 6.31Rjt, and 6.05Rjt for the helical surface and first three helical body modes, 
respectively, at the precession frequency. Thus, we tentatively identify the dominant pattern with 
the first helical body wave and the "hot spots" (at z ~ 20 & 35 Rjt in Figure 11) in the helical 
pattern as caused by wave-wave interaction between the helical surface wave and the first helical 
body wave as illustrated qualitatively by the lower panel in Figure 11. Additional information 
concerning wave-wave interactions is provided by profiles of pressure, axial velocity, and transverse 
velocity components along the z-axis shown in Figure 12. These profiles reveal that the transverse 
velocity components undergo out of phase sinusoidal oscillation (indicating helical twisting at jet 
center) which dies away and then returns at z ~ 55Rjt. The difference between the wavenumbers 
of the first and second helical body modes (see Table 3) can lead to a beat pattern with wavelength 
of ~ 120Rjt similar to that revealed in Figure 12. While the pressure and axial velocity remain 
constant out to z = 25Rjt, Figure 12 shows the development of quasi-periodic pressure and axial 
velocity oscillation with an observed wavelength ~ 6.7 Rjt at larger distance that can only be 
produced by the second pinch body mode. 

We show that the hot spots can be fit by interaction between the helical surface (Hs) and 
first helical body (Hbl) waves, and that the structure at jet center can be fit by a beat pattern 
between first (Hbl) and second (Hb2) helical body modes combined with the presence of the 
second pinch body mode (Pb2) in Figure 13. In Figure 13 fluctuations were computed using 
equations 6 - 8 in §3.4. Component structure is shown along ID cuts parallel to the jet axis 
(z-axis) at different locations on the transverse -l-y-axis where velocity components Vx and Vy 
represent azimuthal and radial velocities, respectively, like Figure 10. The pressure and velocity 
fluctuations that would accompany the helical surface (Hs), first body (Hbl) and second body 
(Hb2) waves using solutions computed from the dispersion relation for the Mach 20 DM jet at 
the precession frequency ooR/u = 1 are shown. This frequency is above the resonant frequency for 
both body modes and lies on the high growth rate plateau for the surface mode wave. The axial 
phase shift in pressure and radial velocity fluctuations between jet center and surface for Hbl is 
much smaller than that shown in Figure 10 for the Mach 5 solution. The difference is entirely the 
result of the difference in the frequency relative to the resonant frequency for this body mode. 
The axial phase shift in pressure and velocity fluctuation between jet center and surface for Hb2 
is similar to that shown in Figure 10 for the Mach 5 solution. Note that significant pressure and 
velocity fluctuations associated with the surface wave occur only near to the jet surface (dashed 
line). Thus, the surface wave does not influence the jet interior at this high frequency. This result 
supports the conclusion that at high frequencies the helical surface wave mode can effect only a 



shallow layer near to the jet surface ( Hardee &: Norman 1988 ; Hardee fc Stone 1997 ), whereas the 



body waves can still effect a larger volume of the jet beam. 

Figure 13 also shows the result of wave-wave interaction between the surface and first two 
helical body modes (Hsbl2). The beat wavelength (best illustrated by the transverse velocity 
fluctuations along the jet axis) is comparable to that seen in the simulation. Additionally, the 
figure shows the pressure fluctuation at the jet surface (dashed line) associated with the wave-wave 
interaction between the surface and first two helical body modes. Here we note an additional 
beat wavelength that arises between the surface and body modes near the jet surface, where the 
pressure fluctuations achieve maxima at axial distances z/R ~ 20, 42 and 65. This confirms 
the tentative identification of the hot spots in the simulation with interaction between surface 
and body helical waves. Note that the beat wavelength shows a slow pressure variation along 
the jet axis (dotted line). A rise in the pressure along the axis occurs as the transverse {vx,Vy) 
velocity fluctuations decrease. This behavior is identical to that seen in the body mode wave-wave 
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interaction (Hbl2) shown in Figure 10. There is a small accompanying axial velocity variation on 
the axis. Just off the axis a short wavelength pressure and axial velocity fluctuation (solid line) 
is evident. This short wavelength fluctuation can couple to the pressure and velocity fluctuation 
expected to accompany the second pinch body mode (Pb2) to produce the short wavelength 
quasi-periodic pressure and axial velocity oscillation seen in the simulation with an observed 
wavelength ~ 6.7Rjt- Thus the simulation results again show that a wave-wave interaction 
between helical body modes can excite a pinch mode. 



5. Summary 

5.1. Comparison with theory 

We have found that the analytical linear analysis provides both good qualitative and 
quantitative agreement with our 3D simulations. The linear analysis successfully indicates: (1) 
the dominate modes in various perturbation frequency ranges, (2) the increase in instability due 
to cooling, and (3) the frequency range in which mode-mode coupling is most likely to occur. 
Quantitative measurements of parameters such as the wavelengths of individual modes and the 
beat patterns formed by wave-wave interaction are correctly predicted and the order of occurrence 
of higher order surface modes are also matched quantitatively with the predictions from the linear 
analysis. In addition, velocity and pressure fluctuations computed from the theory appear to be 
in good quantitative agreement with fluctuations observed in the numerical simulations. The 
coupling between normal wave modes that we have found here has also been found in numerical 
simulations using a spectral type code by Keppens & Toth (1999) who also have identified 
analytically the quasi-linear coupling mechanism that connects a driven helical mode to several 
other low order normal modes. 



5.2. Disruption of the jet beam 



3D simulations confirm results and conclusions reached in 2D simulations (Stone, Xu, & 
Pardee 1997 ). Cooling (especially DM type cooling) increases the linear growth rate of unstable 



modes, especially the helical surface and body modes, and helps to break up the jet beam 
particularly at perturbation frequencies at or below resonance relative to the adiabatic jet. 
What is unique in 3D is that higher order K-H modes play important roles in the jet evolution, 
filamentation and disruption. Animations of our 3D simulations and Figures 2, 5, 6 and 7 show 
that disruption of the jet beam can be categorized by two different fundamental processes: (1) 
the amplitude of distortion caused by helical and other low order surface wave modes grows large 
and the jet beam breaks into filaments in the direction of jet flow. The increase in the contact 
area between the jet beam and the ambient gas leads to enhanced mixing and leading edge shock 
and knot formation, and (2) the compression wave front established inside the jet by the initial 
perturbation is not reflected back into the jet beam from the jet-ambient interface but deposits 
energy at the jet surface and opens a gap in the jet surface, e.g., by a damped helical body mode. 
The gap creates a shock in the ambient gas, and the jet is disrupted. 

It is clear that development of jet distortion, knot formation and jet disruption depend on 
the radiative cooling rate and the perturbation frequency. When the development of jet distortion 
is slowed, by whatever process, disruption is delayed. In general, radiative cooling speeds those 
processes leading to jet disruption relative to an adiabatic jet. DM cooling with a relatively 
shallow dependence of radiative cooling on temperature in both jet and external fluid is somewhat 
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more destabilizing than MB cooling with a steep dependence of radiative cooling on temperature 
in both jet and external fluid. Jets perturbed at high frequency even though helical surface 
and body modes have large growth rates are relatively stable to disruption. It is the value of 
the perturbation frequency relative to the resonant frequency and not the absolute value that is 
important. Typically, our jets remained collimated to twice the distance when perturbed at a high 
frequency when compared to jets perturbed at the resonant or lower perturbation frequency. 



5.3. Knots & Shock Spurs 

We note that knots in the present 3D simulations were less dense than in comparable 2D 
simulations. Mode-mode interactions, also present in 2D simulations, are more complicated in 3D. 
The interactions, such as helical surface and helical body mode interaction and coupling to pinch 
modes, lead to new mechanisms for the formation of staggered emission knots along the jet beam 
which can lead to staggered bow shocks or lead to relatively stationary knots in the jet beam. 

(1) Shocks outside the jet beam created by large amplitude surface displacement provide one 
mechanism for knot formation and for the formation of shock spurs. In this case, the jet-ambient 
interface creates a spiral shaped shock around the jet beam in the ambient gas. If viewed edge-on, 
this shock presents a staggered pattern of shock spurs down the jet beam. Interior to the jet 
surface a strong density enhancement near the surface of the jet beam provided by large amplitude 
helical surface displacement accompanies the spiral shaped shock. Viewed edge-on the density 
enhancement looks like localized knots at the base of the shock spurs. These knots and shock 
spurs move with the speed of the helical surface wave. 

(2) Wave-wave interactions provide a second mechanism for knot formation. In general, our 
results indicate some suppression of knot formation at low precession frequencies resulting from 
the absence of multiple body waves and the associated wave-wave interactions. In wave-wave 
interactions, the knots are formed internal to the jet beam and in one simulation the knots are 
seen to be stationary in temporal animations. Thus, no knot driven shocks are created in the 
ambient gas. Note that the material moving through the stationary knots is moving at nearly 
the jet speed and large radial velocities could still be observed from the spectral lines. However, 
in another simulation the knots formed as a result of wave-wave interaction near to the inlet are 
initially nearly stationary, but at larger distance the knots move and develop shock spurs. At high 
frequencies surface distortion leading to shock spurs does not occur. 

(3) As a jet breaks up mixing with the external environment leads to shocks and knot formation. 
Knots formed in this fashion are primarily associated with the bow-shock region at the head of a 
protostellar jet. In our simulations this region moves outwards into an already outwards moving 
region where jet material previously existed. As a result the resulting bow shock is weaker than 
might be the case otherwise. 



6. Discussion 



Internal knots are one of the most prominent features of YSO jet s. Knots have bee n seen 
to form at the source a nd move out wards with time, e.g., HH 80/81 ( [Marti et al. 1995 ) and 
the Serpens Radio Jets ( Curiel 1995| ). Since they bear many of the same spectral signatures as 
individual HH objects, they should be related to radiative shocks. In our study, all the initial 
perturbations are chosen to be fairly small so as not to excite nonlinear dynamic effects. As a 
result, while we see complex internal structures we do not see strong shocks internal to the jet 
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beam as long as the jet is not disrupted. Biihrke, Mundt, & Ray (1988) and Bodo et al. (1994) 
have proposed that the K-H pinch instabihty may cause knots to form in the jet beam. Our 
simulations do provide evidence for knots produced by pinching at wavelengths induced by helical 
wave-wave interaction. In our simulations the knots formed at the beat wavelength typically 



have negligible proper motion. There are pulsating jet models (see Stone &: Norman 1993b and 
references therein) in which strong internal shocks are created in the jet beam. Knots created by 
the pulsating jet will have proper motions that are comparable to the jet speed. 

Evidence for small and/or accelerating knot proper motion is provided by detailed emission 
line observations of HH 83 by Reipurth (1989). In this highly collimated jet there are about ten 
emission knots, and the jet terminates at a strong bow shock. The distance to the source is 460 pc 
(the source is in the Orion Cloud), and the inclination angle of the jet is ~ 45°. Thus, the physical 
extent of the jet is about 0.1 pc. From the observed width of the brightest knots, the physical 
length-to-width ratio is at least 17. Therefore, from the source to the terminating bow shock the 
linear scale is at least SORjt- Detailed spectral studies reveal that there is a systematic increase of 
the velocity from knot A to knots I and J from —80 km s~^ to —180 km s~^, and the bow shock 
has Ha line emission only, which indicates a shock with velocity of 80 km s^^ or greater. 

These observational results may be compared to our numerical simulation of a Mach 5 DM 
cooling jet (jet speed of about 60 km s~^) perturbed with a precession frequency uRjt/u = 0.1 (a 
period of about 840 years) shown in Figure 8. This volumetric rendering of p^T at an angle of 
inclination of 45° is at simulation time Tsym = 100 (a time of about 6700 years). In the simulation 
the linear scale of the knot forming region is < lOORjt (< 0.1 pc). In the simulation the "hot 
spots" occur at the intersections of helical surface and body waves. The intersections interact 
with the external medium and are shocked as the jet moves supersonically through the ambient 
gas, and shock spurs develop in Figure 8 at knots farther from the origin. Hot spots closer to the 
nozzle are relatively "younger" and are just formed. When a hot spot is newly formed as a result 
of wave-wave interaction it should be nearly stationary. However, the hot spot material moves 
with the jet flow and generates a strong shock at the jet-ambient interface (the shock spurs in 
Figure 8) the material behind the shock spur should have a velocity of about 1/4 of the shock 
speed (also about 1/4 the jet speed). Observed emission should come from the higher temperature 
shocked material at the jet ambient interface. As time passes, the knot is accelerated by the jet 
and the gas in the knot mixes with the jet gas. Thus, an observation should show that a knot 
which is closer to the nozzle has lower velocity, and knots farther from the nozzle, i.e., accelerated 
by the jet, should move faster. Note also that some "hot spots" are on the far side the jet beam, 
which, as the helical wave develops, move away from us along the line of sight with respect to the 
jet beam while other hot spots are on the near side of the jet beam and move toward us along the 
line of sight relative to the jet beam. This mechanism provides additional cyclic variation in the 
radial velocity of knots. 

The initial transverse velocity perturbations used in our simulation are quite small and cannot 
account for the large velocity variation among the knots observed in HH 83. However, the effect 
would be larger if the initial perturbation had a considerable amplitude. The termination of the 
observed jet after knot J would result from disruption of the jet beam. Further downstream, the 
interaction of the jet material with ambient gas would form the classic bow shock - structure 
associated with the working surface observed in many HH objects. Because the morphology of 
the source and of the simulation are very similar we believe that we have found a mechanism that 
could result in variation or systematic increase in radial velocity of knots. 

Another example of similarity between observational morphology and our simulations is 
provided by the structure of the jet in HH 111. In particular, the [S II] images of the jet reveal 
a spiral emission feature wrapped around the jet beam with embedded knots, and the Hq, images 
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present us with remarkable staggered bow shocks farther out along the jet beam (Reipurth et al. 



1997). We can immediately identify the spiral emission feature and bright knots with features like 
those seen in the Mach 5 DM cooling jet perturbed with a precession frequency of toRjt/u = 0.4 
(a precession period of about 210 years) or with the Mach 20 DM cooling jet (jet speed of 236 km 
s~^) perturbed with a precession frequency of LoRjt/u = 1 (a precession period of about 21 years). 
Volumetric rendering of p^T of the Mach 5 and Mach 20 simulations are shown in Figures 5 & 11, 
respectively. In particular, helical emission bands and the knots seen in the Mach 20 simulation 
and in the [S II] map from Reipurth et al. (1997) have wavelength and spacing relative to the jet 
radius that is similar in both the simulation and the observed jet. That this type of structure 
appears in both simulations at very different precession periods and jet velocities suggests that no 
fine tuning is required to generate such structure. Because of the limited domain size in the Mach 
20 simulation, we do not know whether the observed staggered bow shock can be formed at larger 
distance. However, the Mach 5 simulation that develops faster spatially suggests that this would 
be so. 
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Table 1: Physical and numerical parameters used in 3D simulations. 



Mach # Type ^ Grid 7./Rjt 




ioRjt/u 6{rad) XYzones/ Rjt Tsim 


20 A,D 200 X 80 X 80 500 


±4 


0.025 0.0025 12 


25 


20 D 100 X 80 X 80 60 


±4 


1.0 0.01 10 


5 


5 D 200 X 94 X 94 400 


±8 


0.1 0.0075 10 


iUU 


5 A,D,M 200 X 80 X 80 60 


±4 


0.4 0.03 10 


1 

10 


5 A,D,M 200 X 120 x 120 300 


±10 


1.0 0.075 6 




^ Cooling type - A: Adiabatic, D: DM cooling, and M: MB cooling 




Table 2: Growth rates of surface modes for Mach 5 jets precessed at ooRjt/u = 0.4 




Mode kiRjt (AD) kiRjt (DM) kiRjt (MB) 


2TT/ki (AD) 2TT/ki (DM) 


2iT/ki (MB) 


Helical ^ -0.062 -0.087 


-0.068 


imRjt 72Rjt 


92Rjt 


Elliptical 2 -0.135 -0.273 


-0.131 


46Rjt 2mjt 


48Rjt 


Triangular ^ -0.229 -0.274 


-0.210 


27Rjt 23Rjt 


30Rjt 


Rectangular ^ -0.336 -0.327 


-0.303 


ISRjt 19 Rjt 


2lRjt 


^ Growth rate at precession frequency ooRjt/u = 0.4 
^ Maximum growth rate 






Table 3: Helical S, Bl & B2 modes on jets precessed at uRjt/u = 1.0 




Type knRjt (S) knRjt (Bl) knRjt (B2) 


AkRR/t 


kiRjt (S) kiRjt (Bl) kiRjt (B2) 


Mach 5 


Jets 






AD — 1.30 1.46 


0.16 


— -0.072 


-0.054 


DM — 1.18 1.42 


0.24 


— -0.043 


-0.025 


MB — 1.30 1.44 


0.14 


— -0.080 


-0.032 


Mach 20 


Jet 






DM 1.3059 1.0339 1.0868 


0.0529 


-0.0572 -0.0012 


-0.0035 



1 [AkRRjt = kRRjt (B2) - kRRjt (Bl)] 



7. Figures 



10" 
10" 

10" 
10" 
10" 

10" 

10" 
10" 
10" 



— I — I I I 1 1 ii| 1 — I I I 1 1 ii| 1 — r-TTTTFi 

Surface Modes (Ad),--'' 



Rs 



Hs 



I I I lllll | 1 I I lllll | 1 I I I JJjj 

Pinch (Ad) 



.B3 - . - 

-Bl'' 
' S 



-HH 



h\ 1 I I lllll| 1 I I I I II! 



Helical (Ad) 




<jR/u 




— I — I I I I iiij 1 — I I I I iiij 1 — r—rprri 

Surface Modes (DM,).''' 



I I I m il — I I I m il — I I I l iii-t 



Pinch (DM) 



n I m i if 1 1 I I iiiii| — I \ \ nui 

Helical (DM) J^'"^ 

B3 ;;; = -;>-''" 



iTl B2 B3 



I I I lllll f' I I I lllll| — I I l ^ujlt 
Elliptical (DM) 



s ->,--•- 



_B3- 
r-Bi' 



B2 B3 



Hi 



I I I I I ivi 1*^ 



10" 



10" 



10" 



10" 



— I — I 1 — I 1 — r-r-rrpi 

Surface Modes (MB),,''' 



Rs 



Hs 



I I I M ill 1 I ' lll l ll 1 I I 1 1 l li 



Finch (MB) 



-BI " 
■ S 



n I m ii p MM m ii| ' i i i i jji^ 

Helical (MB) 

' S > 



Elliptical (MB) 



-f-H+ttt 



. -B3 - 
■-HI ' 



III Lll I I I I ill if LJ 



10" 



10" 



10" 



cjR/u 



cjR/u 



Fig. 1 — Surface wave solutions to the dispersion relation for the helical (Hs), elliptical, triangular 
and rectangular (Rs) normal modes on the Mach 20 Ad(iabatic), DM cooling and MB cooling 
jets (top panels). Surface (S) wave solutions and first three body wave (Bl, B2, B3) solutions 
to the dispersion relation for the pinch, helical and elliptical normal modes on the Mach 20 
Ad(iabatic), DM cooling and MB cooling jets (lower panels). The dotted lines give the real part 
of the wavenumber, fc^, and the dashed lines give the absolute value of the imaginary part of 
the wavenumber, as a function of the angular frequency, u. Quantities are scaled by the jet 
radius, Rjf, and the jet velocity, u. 
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Fig. 2 — Volumetric rendering of density for the Mach 20 DM jet (top) and Adiabatic jet (bottom) 
with ujRjt/u = 0.025. The volume shown is 500i?jt x ^Rjt x ^Rjt- The images are underscaled 
horizontally by a factor six. Darker shading indicates higher density. 




Fig. 3 — Entrainment volume (top) and entrained mass M§ (bottom) vs. axial position for a 
Mach 5 DM jet with uRjt/u = 0.1. Lines are at values of e = 0.01, 0.02, 0.05, 0.1 & 0.5. and 6 = 
0.5, 1, 2 & 5. 
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Fig. 4 — Volumetric rendering of temperature (top) and density (bottom) for a Mach 5 DM jet 
at ujRjt/u = 0.4. The maximum temperature and density are ~ 40,000 K and ~ 1,400 cm~^, 
respectively, and the volume shown is Q^Rjt x ^Rjt x ^Rjt- Darker shading indicates higher 
temperature. In general higher densities are darker but in the first half of the jet lighter areas are 
also above the initial jet density. 




Fig. 5 — Composite p^T volumetric rendering of the Mach 5 DM jet at ujRjt/u = 0.4 shown in 
Figure 4. The jet is inclined at 45" to the line of sight and the volume shown is 60i?jt x SRjt x ^Rjf 
Darker shading indicates higher values of p^T. 



Fig. 6 — Center cuts in density of Macli 5 AD (bottom), DM (middle), and MB (top) jets at 
ujRjt/u = 1.0. Mapping to greyscale is the same in each panel and the panel size is 300Rjt x 20Rjt 
and underscaled horizontally by a factor of 3. In general, lighter shading indicates higher density 
but in the jet the darkest shading indicates the highest densities. Densities range from ^ 60 cm~^ 
in the ambient to a maximum of ~ 1, 200 cm~^ in the jet. 
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Fig. 7 — Mach 5 DM jet cross sections at uRjt/u = 0.4. The number in each frame indicates the 
distance from the inlet in units of R^. In general, a lighter shade in these cross sections indicates 
a higher density but in the jet the darkest shading within lighter regions indicates the highest 
densities. Densities range from ~ 60 cm~^ in the ambient to a maximum of ~ 1,200 cm~^ in the 
jet. 
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Fig. 8 — Composite p^T volumetric rendering of the Mach 5 DM jet at coRjt/u = 0.1. The jet 
is inclined at 45° to the line of sight. Panels show the entire jet (top), and entire jet excluding 
material with T < 1000 K (middle) with volume 300Rjt x 16Rjt x IQRjf. Lower panel shows the 
knot region excluding material with T < 1000 K with volume bORjt x IGRjt x IGRjt. Darker 
shading indicates higher values of p^T and reveals the hot dense parts of the jet. 
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Fig. 9 — Plots of pressure (top), axial velocity (middle) and transverse velocity components 
(bottom) along the jet axis from the Mach 5 DM jet at uRjt/u = 1.0. 
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Fig. 10 — (top two rows) Pressure, axial velocity (v^), azimuthal velocity (v^) and radial velocity 
(fy) for the first two helical body modes Hbl and Hb2 on a Mach 5 DM jet along ID cuts parallel 
to the z-axis at radial locations rjR = 0/8, 1/8, 2/8, 7/8 on the +y-axis. The outermost 
(innermost) radial locations are indicated by the dashed (dotted) lines, (bottom two rows) 
Pressure, axial velocity (vz) and transverse velocity components for wave- wave interaction between 
the two helical body modes (Hbl2) along the jet axis, and for the pinch cooling mode (Ps2). For 
Ps2, Vr is shown on the axis (dotted line) and near to the jet surface (dashed line). 
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Fig. 11 — Composite p^T volumetric rendering of the Mach 20 DM jet at ooRjt/u = 1.0 (top) 
and illustration of the interaction between surface and body wave (bottom). The size shown is 
GORjt X SRjt- Darker shading in the volumetric rendering indicates higher values of p^T. 
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Fig. 12 — Plots of pressure (top), axial velocity (middle) and transverse velocity components 
(bottom) along the jet axis from the Mach 20 DM jet at ojRjt/u = 1.0. 
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Fig. 13 — Similar to Figure 10 but for helical surface (Hs), helical body modes Hbl and Hb2 and 
pinch body mode Pb2 on a Mach 20 DM jet. The pressure panel for the combined helical surface 
and body waves (Hsbl2) shows ID cuts along the jet axis (dotted line), at r/R =1/8 (solid line) 
and at r/R = 7/8 (dashed line). The axial velocity panel shows ID cuts along the jet axis (dotted 
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line), at r/R = 1/8 (solid line). Pressure, axial velocity (v^) and radial velocity for the pinch body 
mode (Pb2) are plotted along the jet axis in the lower panels. 



